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Abstract 

A molten carbonate fuel cell (MCFC) is an electro-chemical energy conversion technology that runs 
on natural gas and employs a molten salt electrolyte. In order to keep the electrolyte in this state, 
the cell must be kept at a temperature above 500° C, eliminating the need for precious metals as the 
catalyst. There has been only a limited amount of research on modelling the transport processes 
inside this device, mainly due to its restricted applicability for mobile applications. 

In this work, three one-dimensional models of a MCFC cathode are presented based on differ- 
ent types of diffusion and convection. Comparisons between models are performed so as to assess 
their validity. Regarding ion transport, it is shown that there exists a limiting case for ion mi- 
gration across the cathode that depends on the conductivity for the liquid potential. Finally, an 
optimization of the diffusivity across the cathode is carried out in an attempt to increase the cell 
performance and its longevity. 
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Nomenclature 



b Bruggeman correction, 1.5 

c Concentration, mol/m^ 

Cj Total concentration, mol/m^ 

D Diffusivity (O2 in air), m^/s 2.5 x 10^ 

F Faraday's constant, C/mol 96487 

io Exchange current density, A/m^ 1 x 10^^ 

L Thickness of cathode, m 8 x 10^^ 

n Number of electrons in cathode reaction, 4 

P Pressure 

R Ideal gas constant, J/Kmol 8.314 

T Temperature, K 900 

u Fluid velocity, m/s^ 

a Transfer coefficient, 0.5 

eg Gas porosity, 0.4 

€/ Liquid porosity, 0.3 

Cg Solid porosity, 0.3 

r] Polarization coefficient, V 

K Permeability, m^ 1.9 x 10~ 

fj. Viscosity, kg/ms 2.25 x 10 

v Stoichiometric coefficient 

ai Liquid conductivity, S/m 140 

as Solid conductivity, S/m 1300 

(p Potential, V 



1 Introduction 



Molten Carbonate Fuel Cells (MCFCs) are widely being considered for stationary power genera- 
tion and a better understanding of the transport processes in the electrodes and cells are needed 
to improve viability. Recently, a MCFC system was built at Enbridge headquarters in Toronto, 
Ontario. Their system features four separate MCFC stacks where the primary fuel, natural gas, is 
pumped into the system from a pressure let-down station. The excess heat released by the fuel cell 
system is used to heat the adjacent building and preheat the gas before expansion. 

Molten carbonate fuel cells were initially developed with the intention of operating directly on 
coal. The primary fuel currently in use is either coal-derived gases or more commonly natural gas 
[1]. MCFCs are still under development and have not reached market acceptance as a possible 
primary or secondary source of energy. 

The concept of the MCFC is almost a century old with the first patent awarded in 1916 to 
W.D. Treadwell. It was first conceived in Europe in the 1940s as an attempt to convert coal to 
electricity in carbonate media. An initial demonstration was successfully completed by Broers and 
Ketelaar in the 1950s, with the first pressurized stack appearing in the 1980s. Most of our current 
knowledge stems from work done in the 1970s and 1980s [1 . 

Current development concentrates on base-load utility applications as well as dispersed or dis- 
tributed electric-power generation with heat co-generation. Due to low power densities and long 
start-up times, there is limited potential for mobile applications, although MCFCs might be suitable 
for power-trains for large surface ships and trains. 

This work presents three mass transport models of a MCFC cathode electrode to be compared 
for the same parameters. The second section of this work provides a brief introduction to the 
physical and chemical processes in the cathode electrode. The third section focuses on the three 
models being used as well as numerical results. Section 4 presents an analytical resolution to the 
non-existence of numerical solutions found by altering the liquid conductivity parameter based on 
values in White et al. [2L Section 5 presents an optimization model for the mass transport of a 
single species in the MCFC cathode electrode. Conclusions will be drawn based on numerical and 
analytical results in the fifth chapter. Before we begin our analysis, a brief overview of a MCFC 
cathode is provided and models based on mass transport phenomena shall be described. 

2 Cathode 

Inside the cathode, oxygen and carbon dioxide flow in the same direction across the electrode at 
different rates. The three-phase boundary allows for reactions to occur along the length of the 
domain of the electrode. Across the channel interface, only oxygen and carbon dioxide are able to 
flow while the liquid electrolyte is kept from flowing into the channel causing corrosion. In fact, 
the liquid electrolyte distribution is only controlled by capillary pressure. Therefore, the pore size 
must be maintained very carefully during manufacturing. At the electrode/electrolyte boundary, 
the electrolyte concentration is much larger and fills the pores, which helps to avoid any gas leakage 
into the electrolyte assembly. 

The inlet gas at the cathode is mainly composed of N2, O2 and CO2 but contains trace amounts 
of other molecules since the main source of oxygen comes from air. The effects of nitrogen are 
neglected in this thesis since it does not flow within the cathode. Only a few metals are stable 
as a cathode material due to the extremely corrosive nature of the molten carbonate electrolyte 
and currently nickel oxide is in use. Only semiconducting oxides are feasible from a cost point of 
view |lj. The mean pore size of NiO electrodes is about 10 /^m. The smaller pores are filled with 
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Figure 1: Domain schematic for a one-dimensional cathode electrode. 



the electrolyte to form the three-phase boundary needed for the reaction, while the larger pores 
remain open for gas flow. Nickel oxide is also slightly soluble in the electrolyte which limits the 
lifetime of the cell, according to 

NiO + C02 ^ Ni2+ + C0j. (1) 

The optimal cathode performance depends upon the gas composition where there exists a 2:1 ratio 
of CO2 to O2 consumed in the overall electro-chemical reaction. In order to reduce the NiO solubility 
and increase lifetime, the carbon dioxide concentration should be reduced, although if it is too low, 
the dissociation of carbonate ions becomes significant 

C0i=^C02 + 0=, (2) 

thereby limiting the cell lifetime due to electrolyte losses. The balance between NiO solubility and 
dissociation of COl" can become very difficult to control, and to predict cell lifetime is generally 
challenging. 

The optimal thickness of the electrode, which depends upon the gas composition and current 
density as well as other operating conditions, ranges from 0.4 — 0.8 mm |T]. 

At the three-phase boundary in the cathode, oxygen and carbon dioxide diffuse towards the 
electrolyte which has penetrated the NiO pore. Where the gas flow meets the electrolyte, the gas 
molecules are absorbed into the electrolyte, react with the electrons at the surface of the electrode, 
and produce carbonate ions. The electro-chemical reaction is given by 

^02 + C02 + 2e- ^ CO^. (3) 
3 Mathematical Models of Diffusion 

The diffusion of gases across the cathode of a MCFC are now studied mathematically using three 
different models. The first model considers only diffusion, while the second and third models 
consider two different types of diffusion as well as convection. 

The electro-chemical reaction in the cathode. Equation ([s]), involves a reaction between three 
constituents and the production of another. The reactants diffuse across the cathode from the 
channel to the electrode/electrolyte boundary (left to right across the domain shown in Figure [T|, 
while the ions also move from left to right. Therefore, a system of four equations will be used that 
describe the flow of gases, electrons, and ions. 

3.1 Fickian Diffusion 

The first model is derived using Tick's Laws of Diffusion and Ohm's Law. 
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The solid and liquid potentials are given by the change in current density using Ohm's Law and 
the Butler- Volmer equation 



) = -i'sS{(l)s,4'l,Co2,Cco2), (4) 



d 

dx \ " ' dx 

^ ^crie/(x)''^^ = viS{(t)s,(t)uCo^,Cco^), (5) 

where ai and as are the conductivities for the liquid and solid phase, respectively. The solid (e~) 
conductivity in the MCFC is generally one to two orders of magnitude larger than the liquid (COl") 
conductivity, which will create a near constant potential for the former across the domain. The 
potential is described by (j)i and (jjg, and determined by the porosity ei{x) and es(a;), stoichiometric 
coefficients Us and f;, and the Butler- Volmer equation, S, describing the reaction kinetics. 

Using Tick's Second Law, the rate of change of concentration as a result of diffusion by O2 and 
CO2 and their reaction is found to be 

-L>eg(x)^^^ ) = -i/o25'((/)s,</'i,Co2,Cco2), (6) 



dx \ dx ^ 

^ ("^^9(2;)^-^^ = -I^CO2 5'((As,0«,Co2,CcO2), (7) 

where D represents the diffusion coefficient, which is the same value for both gas species. Note the 
exponential dependency on the porosity eg{x). The constant h represents the Bruggeman correction 
for the impedance of diffusion in each phase [3] . 
The Butler- Volmer equation is given by, 

rjanF 

S'(0s,</'i,Co2,Cco2) = «0Co2Cco2e «^ , r] = (j)s - (1)1, (8) 

where r] represents the difference between solid and liquid potentials, n represents the number of 
electrons in the reaction, and io the exchange current density. The negative exponential term in Q 
is dropped since the polarization is at least t] > 0.05 V. 

The correction factor for the diffusivities and conductivities are given by the porosity, e^. The 
volume fractions must add up to one. If the porosity of two quantities is known, the third quantity 
is given by, 

eiix) + eg{x) + es{x) = 1. (9) 

Equations can be solved numerically using Newton's Method for computational efficiency 

or iterative methods. 

3.2 Fickian Convection-Diffusion 

Using a convection-diffusion model, the flux becomes the sum of the convective flux and the diffusive 
flux. In this model, the molecular interactions are not considered explicitly. They will be included 



in the third model (Section 3.3). 



The convection-diffusion equations for O2 and CO2 are given by 

^ ^ucoj - Degix)'''^^^ = -Uo2S{(t>s, (t>i, C02, CC02), (10) 

^ (uCc02 ~ ^'^ai^^^^'-^) = -^C02'S'(</'s>«,C02,Cc02), (11) 
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where the fluid flux, u, is described by Darcy's Law for porous media, 



« = -^^^ = -^iir^. (12) 

ax fi ax 

Compared to the Fickian Difl'usion model, there is an extra term and Newton's method is used to 
solve this type of non-linear problem. Note that the two diffusive fluxes do not necessarily add up 
to zero. 



3.3 Multicomponent Convection-Diffusion 

The Maxwell-Stefan equations are used to describe the flux of the diffusing species as well as the 
interactions between different molecules in the fluid flow. The flux for this model is the sum of 
the convective flux and the multicomponent diffusive fluxes. The conservation equations for the 
concentration are given by 



- (^uco^ ~ ^^'''-^''^^"'ci^T^) ^ -^02S{(t)s,(t>l,Co2,Cco2), (13) 
^nccoj - Deg{x)''cT-^'^^ = -J^co2'S'(0^i, <A/, C02, CcoJ- (14) 



This model is similar to the model presented in White et al. [2], with the addition of convection. 

These equations become highly non-linear due to the total concentration term, Cy, which also 
appears in Darcy's law (cj = + Ccoj)- Newton's Method can be used to solve this model but 
it is more computationally expensive than the Fickian convection-diffusion model. Note that here 
the two diffusive fluxes always add up to zero. 



3.4 Boundary Conditions 

The gas enters the electrode at the channel while the liquid electrolyte penetrates the electrode 
pores at the opposite side. The gas is unable to flow into the electrolyte due to the pore fllling and 
the electrode is manufactured in such a way as to avoid the corrosive electrolyte penetration of the 
gas channels. 



3.4.1 At the Channel {x = 0) 

At the channel, the gas is flowing into the electrode. This boundary uses Dirichlet conditions that 
give the value for the concentration of O2 and CO2 as well as the solid potential. The electrolyte is 
not allowed to move from the electrode into the channel and zero flux is enforced using a Neumann 
condition for the liquid potential. The boundary conditions are chosen as 

In reality, the electrolyte distribution in the cathode is non-uniform and e/(x) will be zero for a 
finite range < x < xo- In fact, e/(x) is usually a monotonically increasing function. In essence, 
we are moving the channel from rc = to x = xq in this work since there are no reactions between 
and xq. 
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Figure 2: The solid and liquid potential dis- 
tribution across the electrode using the reference 
parameters as stated in the nomenclature. The 
Fickian convection-diffusion (Conv-Diff) and multi- 
component convection-diffusion (Binary) models 
share the same profile. 



Figure 3: The potential distribution with the liq- 
uid conductivity decreased from the reference value 
of f40S/m to 50S/m, which increases the reaction 
rate. 



3.4.2 At the Cathode/Electrolyte Interface {x = L) 

At the cathode/electrolyte interface, the gas and electrons are not allowed to enter the electrolyte 
and this is reinforced using Neumann conditions. The liquid potential is given by a Dirichlet 
condition. These boundary conditions are chosen to be 

^(L).o. ML)^,,., ^w^o. '^m... a.) 



3.5 Results 

3.5.1 Solid and Liquid Potential 

The potential drop is shown by plotting the solid and liquid phase potentials in Figures [2]|5j 

The potential difference across the domain using the reference values as in the nomenclature 
is plotted in Figure [2| The potential remains relatively constant across the domain and all three 
models remain within one-thousandth of a decimal point in agreement. 

The steady state solution for the potential difference, as the liquid conductivity is decreased, 
is plotted in Figure [sj The liquid conductivity, ai, is decreased from 140 S/m to 50S/m. If the 
conductivity is decreased further, the steady state solution ultimately does not exist and this will 
be examined in Section |4} As the liquid conductivity is decreased, the potential gradient increases, 
hence, the reaction rate increases and since all four equations are coupled, the solid potential 
and species concentrations also change. This increase is attributed to the inverse relationship 
between the rate of change in liquid potential and conductivity. The convection-diffusion and 
multicomponent models share the same profile despite the differences in the formulation. All three 
models have a change within the same order of magnitude. 

The value for the permeability of the MCFC cathode is difficult to obtain but is comparable 
to that of a Proton Exchange Membrane fuel cell (PEMFC) gas diffusion layer, which contains a 
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Figure 4: The potential distribution with the per- 
meability decreased three orders of magnitude ef- 
fectively turning off the convective flux. 



Figure 5: The potential distribution with the ex- 
change current density decreased by two orders of 
magnitude, slowing the reaction rate. The solid po- 
tential varies less than 10"'^ V across the domain. 



similar pore size. The potential difference, when the permeability, k, is decreased by three orders of 
magnitude, causing the convective flux to approach zero, is plotted in Figure|4j With the convective 
flux near zero, the Fickian diffusion and convection-diffusion models are essentially identical and 
this is shown in the graph as they now share the same profile. The enforcement of zero net diffusive 
flux in the multicomponent model keeps the potential change lower than that shown by the other 
models due to the small mass transfer of gas. If the permeability increases, the convective term 
dominates further and will result in similar results as before. 

The exchange current density, io, enters the Butler- Volmer equation and can be found to vary by 
several orders of magnitude depending on the model. The potential difference is plotted in Figure 
[5] with the exchange current density at a value of two orders of magnitude less than the reference 
value, 1.0 X 10^'^A/m^, which is the value used in White et al. [2]. By decreasing the exchange 
current density, the reaction rate for each species decreases and the rate of change decreases as 
well. At such low reaction rates, each model shares the same profile across the domain. 

3.5.2 Concentration 

The change in concentration of O2 and CO2 is plotted in Figures [6]j9j corresponding to the results 
for the potential difference. 

The convection term dominates the diffusion term for O2 across the electrode as it pushes the 
gas towards the electrolyte as in Figure |6] using the reference values in the nomenclature. This can 
be seen due to the increase in O2 concentration as it approaches the cathode/electrolyte boundary. 
With convection dominating, there is less than a 10% change in concentration, while the Fickian 
diffusion model without convection shows a 10-20% change in O2 and CO2 concentration. 

Considering the concentration of O2, it increases across the domain due to convection and it 
is not possible to approximate the convection-diffusion profile using an effective diffusivity for the 
Fickian diffusion model as in White et al. [2J, which did not include convection. 

The steady state solution for a decrease in the liquid conductivity, ai , is shown in Figure [7| The 
value is decreased by an order of magnitude and if decreased further, the steady state solution does 
not exist (see Section |4| and the code no longer converges. Since the system of four differential 
equations is coupled, the gradient of the concentration of the gas-species changes along with the 



8 



Concentration „ Concentration 




Figure 6: The O2 and CO2 concentration pro- 
files across the electrode using the reference val- 
ues. The Fickian convection-diffusion (Conv-Diff) 
and multi-component convection-diffusion (Binary) 
models share the same profile. 



Figure 7: The concentration profiles with the liq- 
uid conductivity decreased from the reference value 
of 140 S/m to 50S/m, which increases the reaction 
rate for the liquid potential. 



changes in reaction rate. The convective flux remains dominant and the convection-diffusion and 
multicomponent models continue to share the same profile. 

The concentration profile is shown in Figure [8] when the permeability, k, is decreased by three 
orders of magnitude. Since the convection is close to zero, diffusion dominates and the Fickian 
diffusion and convection-diffusion models are essentially the same as was described previously for 
the potential difference. The differences between the convection-diffusion and multicomponent 
models becomes more apparent as the concentration drops almost 80% across the cathode using 
the latter model, while less than 20% for the former model. This can be attributed to maintaining 
a non-zero convective flux. The discrepancy between the two models will be further analyzed in 
Section 13.5.41 

The exchange current density, zq, when decreased by two orders of magnitude, decreases the 
reaction rates for all species. As the reaction rate decreases, less species react at the three-phase 
boundary and the concentration drop is minimal across the domain. With the permeability at the 
standard value, convection dominates even with low reaction rates. Once again, this is the value 
of the exchange current density reported in White et al. |2]. It still shows the convection term 
dominating which must be considered in the mass transport of the MCFC cathode. 

It should be pointed out here that the model by White et al. [2] is ill-posed in that it neglects 
convection while keeping binary- diffusion terms. Since the latter add up to zero, we cannot have a 
net flux of gas across the electrode, an obvious contradiction, which they do not discuss. This can 



be seen by adding (13) and (14) giving (for n = 0), 



(^-Deg{x)^CT^ (^ + ^)) = -i''02 + ''c02) S{<l)sAl, Co,, CC02), 



^ec;(^)'cT:^(l)) =0, (17) 



which is a contradiction since S is not identically zero. 
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Figure 8: The concentration profiles with the per- 
meability decreased three orders of magnitude, ef- 
fectively turning off the convective flux. 



Figure 9: The concentration profiles with the ex- 
change current density decreased by two orders of 
magnitude, slowing the reaction rate. 



3.5.3 Cell Performance 

The MCFC operates at current densities in the range of 100 — 200mA/cm^ with a cell potential 



of 0.75 — 0.90 V. Figure 10 shows the polarization curve for the convection-diffusion model at 
several different values of the diffusivity, D. As the diffusivity decreases, the polarization curve 
loses its linear character as the degree of irreversibility increases. At the standard reference value 
for the diffusivity, 2.5 x 10^^ m^/s, the cell operates in the appropriate range for the current density 
and potential difference. This graph was made using an open-circuit or reversible cell potential of 
Er = 1.0 V, which is appropriate for MCFCs The expected cell polarization for the half-cell 
reaction in the cathode is shown in the figure as a box. 

3.5.4 Comparison of Convection-Diffusion Models 

In many cases, researchers in the field of porous media flow are interested in approximating the 
convection-diffusion equation by using only Fickian diffusion with an effective diffusivity that will 
provide a similar profile with convection excluded. In the case when the convection term is included 
along with Fickian diffusion, the concentration of O2 can increase from the channel to the electrolyte 
when convection is in the opposite direction of diffusion. Due to this increase, it is not possible 
to use an effective diffusivity that will provide the same profile since diffusion will only cause the 
concentration to decrease across the domain. This means that convection must be included in the 
mass transport of the MCFC cathode and can not be neglected as in White et al. [2^ . 

Another area of interest in the study of mass transport in fuel cells is the comparison between 



the convection-diffusion model of Section 3.2 and the multicomponent convection-diffusion model of 
Section 3.3 Depending upon the parameters, it is possible that Fick's Law can provide a sufficient 
description of diffusion across a domain. Note that in the first model, the diffusive fluxes do not 
add up to zero whereas in the second model they do. 

The permeability, k, is a measure of the ability of a fluid to move through a porous medium. It 
can be very difficult to measure the permeability of a material for fuel cells as it is directly related 
to the manufacturing of the material. Values for permeability in the literature for MCFCs have not 
been found but we can safely use the values for the permeability in the proton exchange membrane 
fuel cells (PEMFC), which is the primary fuel cell in use for mobile applications. Currently, there 
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Figure 10: Cell potential versus current density: polarization curve for the half cell Fickian convection- 
diffusion model. 



exists only 2-3 manufacturers of MCFCs and values for the permeability are hard to come by. The 
pore size in the cathode of the PEMFC is approximately 10 ^m, which is close to that of the MCFC. 
Realistically, we can determine the permeability of the MCFC to be close to 1.9 x 10^^^ m^ as in 
the PEMFC model of Promislow et al. [T. 

The diffusion coefficients in the convection-diffusion model are the same for any number of 
species considered using Pick's First Law of Diffusion. In the Stefan-Maxwell formulation, only n—\ 
fluxes are independent and for n > 2, the diffusion coefficients consist of a matrix. The diffusion 
coefficients for the species in the cathode of the MCFC are also difficult to find from literature. For 
instance. White et al. |2] state the diffusion coefficients to be 1.16 x 10~^m^/s, which is referenced 
from Cussler [5l |2] , but the diffusion of O2 through air, which is the highest possible diffusivity, is 
only 0.25 x 10~^m^/s. The value reported in Fehribach et al. [6], 4.3 x lO^^m^/s, is much more 
realistic. 

Due to the uncertainty of values for the permeability and diffusivities, a comparison between 
the two models for different values is performed so as to validate the use of the convection-diffusion 
equations as an approximation to the Stefan-Maxwell equations. The convection-diffusion model is 
far less computationally expensive than the more non-linear multicomponent convection-diffusion 
model, which will be beneficial for scaling- up the model to higher dimensions. 



Figure 11 shows the difference (error) between the O2 concentration profiles of the Fickian 



convection-diffusion and multicomponent convection-diffusion models. 

The error is calculated as the L'^-norm of the difference between the O2 concentration profiles 
of the Fickian convection-diffusion (fed) and multi-component convection-diffusion (mcd) models 
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Figure 11: The difference between the concentration profiles for both convection-diffusion and multi- 
component transport models. For the standard values of both diffusivity, D, and permeability, k, the 
error between the solutions is minimal. The error increases for smaller values of the permeability where the 
convection term is closer to zero. The colours represent the error between profiles and is less than 10^^ up 
to permeabilities less than 1 x 10~^^, where the error grows above 1. 



normalized by the area under the O2 convection-diffusion profile, so that 



^ l^fcd „mcd 1 2 J 

1^02 ^€02' ""^ 



1/2 



error = -j^ , (18) 

/o dx 

where Simpson's rule was used to find the area. 

For permeability values above 10~^^ m^, there is a very small difference between the profiles (less 
than 10~^). For smaller values of the permeability, the difference increases significantly. However, 
for realistic values based on the known diffusivities from the MCFC literature and permeabilities 
from the PEMFC literature, the convection-diffusion model is in very good agreement with the 
multicomponent model. 

Based on these results, it seems justified to approximate the multicomponent model by using 
the simpler, and computationally less expensive, convection-diffusion model in MCFC cathode 
electrodes. This is beneficial for scaling-up the model to higher dimensions as well as in stack 
models where run time increases dramatically. 
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4 Fickian Diffusion: Non-Existence of Steady-State Solutions 



While performing numerical simulations of the models given in Section 3.1, the solution did not 
converge for small values of the liquid conductivity. The value for liquid conductivity, ai , reported 
in White et al. [2] is around 2.0 S/m and in Fehribach et al. [6] it is around 140 S/m. The model 
presented here will only admit solutions when the conductivity is well above the value from White 
et al. 

The question that arises is whether this is a numerical issue or an ill-posedness of the model. 
Using Equation ([5]), 



d ( .bd-4>l\ ■ f naF{(ps - (pi) 



(^i^iK^) -JZ ) = ^^i^oCojCcoj exp ( — ), (19) 



we can solve this analytically by non-dimensionalizing the equation according to 

4>i,xx - 5e-^' = 0, (20) 



where 



naF - _ r-i X _ ""^03=-^^^ RT 



-Ti- jjrr-Pu x = L X, 6= y-^ ^^0Co2Cco2e'f'^ (21) 

RT aiei[x)° naF 

In this non-dimensionalized model, the porosity ei{x) and solid potential (ps are considered to be 
constant across the domain, which is a valid assumption based on the results of Chapter [3} The 
concentration of O2 and CO2 are also considered constant across the domain which is a reasonable 
first-order approximation. 

The parameter 6 contains the liquid conductivity ai. We will now show that a critical value, 
6c, can be found beyond which the steady-state solution will not exist. 

4.1 Thermal Runaway 



The potential equation, Eq. (20 1, is similar (replace (pi — > to the steady-state equation for 



thermal runaway found in Fowler [7] 

+ Ae^ = 0, (22) 
with boundary conditions 

^(±1) = 0. (23) 
This system can be solved analytically as follows: multiply both sides by 0^ and integrate, yielding 



j O^e^^dx = j -Xe^e^dx, (24a) 

^ ]^el = -Xe^ + C. (24b) 

The constant C can be found by imposing the symmetry condition 9{—x) = 6{x) and if 9 is smooth 
enough, 6*^.(0) = 0. We find C = Ae^'" where 6*0 = 61(0). 
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Considering only the problem on the interval x = [0, 1] in comparison with Eq. (20), the ODE 
is re-arranged and the method of separation of variables is applied which gives 



de 



/6»o V e^o — e" Jo 
Let us set z = V e^o — so that 
dd ■- 



2A\/e^o - e^, 

2Xdx = V2Xx. 



de 



Ize-^dz, 
-2dz 



Substitution into (25b) gives 

r- /"^ de 
J 00 V e^o - 

where we have used the integral |H] 

du 1 ^ -\ u 
— tanh — . 
a a 



a? — 



Re-arranging for z leads to 



dz 



-2e"^o/^ tanh" 



(25a) 
(25b) 

(26a) 
(26b) 



(27) 



(28) 



(29) 



using the identity tanh(— x) = — tanh(x). Defining 7 = \/ ^e^ox, we obtain 



or 



= e^" (1 - tanh^ 7) = e^" sech^ 



e = e^ — 2 In cosh 7, 



7 



(30a) 



(30b) 



using the identity 1 — tanh^ 7 = sech^ 7 = cosh^^ 7. This solution for e was found in Fowler [7. and 
we can use the same approach to solve for the liquid potential. 

For this model based on thermal runaway, the maximum temperature, occurs at a; = (by 
symmetry), and is determined by satisfying the boundary condition at x = 1, 0(1) = 0, and so 



from (30b) 



= cosh 1 \j -e 



9o/2 



(31) 



The solutions to this transcendental equation are studied in Figure 12 which was re-created based 
on the figure presented in Fowler [7]. Based upon these results, there will either be 2, 1 or 
solutions depending on whether A < Ac, A = Ac or A > Ac, respectively. The critical value, Ac, is 
found in the following manner. 
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Figure 12: In order to obtain solutions to Eq. (31 ), the possible solutions based on different values of A are 



shown graphically. As A decreases, the number of solutions increase [7]. 
Define 



1/2 



It; = I - I e 



as well as the functions 
y\ = coshtD, 



y2 = e 



1/2 



w. 



(32) 

(33a) 
(33b) 



There exists only one solution if ?/2 is tangent to yi at some point w = w* . The tangent is given 
by 



d (2V'^ d 
— \ — ] w = — cosh w, 
dw \ X dw 



1/2 



sinh tf, 



(34a) 
(34b) 



and point w* is found from 



1/2 

coshw* = , 
2 V 2 



1/2 



1 + 



1/2 /;, + 2^1/2 



(35) 
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Figure 13: Response diagram for Oq as a function of A. As the value of A decreases, the value of Oq has 
multiple possibilities that will lead to a steady state solution based on the outcomes of Figure [12] [7]. 



Upon substitution of (35) into (34b), the critical value, Ac ~ 0.878, is given by the unique value 
that satisfies 



1 



sinh 




which was found using Maple. 



Figure 13 gives the value of ^-s a function of A. There exists an asymptote at Ac 
are 2 solutions for A < Ac, 1 solution for A = Ac or solutions for A > Ac. 



(36) 



where there 



4.2 Fickian Diffusion 



In order to find a solution to expression (20), let d 



solved analytically using the same approach as the one above. 

This time, applying the boundary condition at the right boundary {(p{l 
value (po which satisfies 



/ in Eq. (22). Then this system can be 
0) gives the minimum 



cosh 




-4>o/2 



(37) 



Once again, we can find the critical value, 6c ~ 0.878, by examining the solutions to Eq. (37) in 
the same way as Eq. (31). Here, 5c is given by 



1 



■ sinh 



(38) 
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Figure 14: Response diagram for (/>o as a function of S. The result is the same as for thermal runaway where 
steady state solutions appear when S < Sc- 



The solutions are the same and given in Figure 12 where 9 = —cj) and X = 6. 

Figure 14 shows the value of 0o as a function of 5. There exists an asymptote at 6c, where there 
are two solutions for 6 < 5c, one solution for 6 = 6c or no solutions for 6 > 6c- In the case of two 
solutions, the "high- volt age" solution (i.e., larger value of (po for given 6) is very likely unstable 
(without proof). Therefore, this solution cannot be observed experimentally. 

In the MCFC cathode model, as the liquid conductivity decreases, 6 will increase (see Eq. (21 )). 
The analytical results presented here explain the non-existence of steady state solutions found in 
the numerical computations. For large values of D and as, resulting in near constant functions 
C02) Ccoj a-iid 4>s, we find very good agreement between the theoretical value of 6c and the value 
obtained numerically. 

Also, as 6 ^ 6c, the CPU time increases dramatically. 



5 Optimization of a MCFC Cathode 
5.1 Full Model Optimization 

Initially, the optimization of the porosity in the cathode was attempted for the full model of four 
variables by using the Fickian Diffusion model with the built-in optimization routines in MATLAB 
(i.e., fmincon). The gas porosity was taken to be a function of the position across the cathode, which 
affects the results through the Bruggeman term (ej(x)''), while the liquid porosity was updated at 
each position based upon a constant solid porosity. The porosity was updated after convergence to 
a solution based upon maximizing the current density at the channel. The flux at the channel for 
the solid potential was used for the maximization criteria. 
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Results show that either the model is ill-posed or there exists a point where the optimization 
routine can no longer alter the porosity to reach a maximum for the current. Each time the 
program was run for the same parameters but different, random initial porosities, the results 
were inconsistent with each other and were not representative of what was expected based on 
the physics of the problem. At the channel, it is expected to have a larger porosity than at the 
electrode/electrolyte interface for the gas porosity to allow for more of the gas species to flow 
into the electrode and react accordingly. Hence, this optimization approach using fmincon was 
not pursued any further and a different optimization problem was investigated instead, employing 
optimal control methods. 



5.2 Optimal Control 

Optimal control problems can be used to measure how effective a given control of a system is by 
minimizing a cost functional. This type of problem is an important part of optimization and has 
many applications. The following is a formulation of optimal control problems from Pedregal |9j. 
The state of a given system is described by a number of parameters, 

X = {xi,X2, ■ ■ . ,Xn), (39) 
which evolve according to a state equation, 

x'{t) = f{t,x{t),u{t)), (40) 
with boundary or initial conditions 

x(0) = xo, x{T) = XT. (41) 

The control parameters are given by 

U = {ui,U2,...,Un), (42) 

and they may depend on t also. 

The functional we are attempting to minimize is defined as 

I{x,u)= [ F{t,x{t),u{t))dt. (43) 



Jo 

Both the state equation and the objective functional depend upon the control parameters u. 
A pair {x,u) is said to be feasible and admissible if the following is fulfilled [9]: 

1. constraints on the control: u(t) E K for all t € (0,T), where K is the permitted range; 

2. state law: x'{t) = f{t,x{t),u{t)) for all t £ (0,r); 

3. end-point conditions: x(0) = xq, x(T) = xt- 

An optimal control problem may consist of having both or only single endpoint conditions, and 
transversality conditions may be needed to complete the formulation of the problem, which will be 
discussed later. 

An admissible pair (X, U) is sought such that, 

I{X, U) < I{x, u) (44) 
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for all other feasible pairs (x,n). 

We can incorporate the point-wise constraint Eq. (40) using a Lagrange multiplier or co-state 
p{t) and consider the augmented functional 

I{x,u,p,x')= [ [Fit,xit),uit))+p{t)-{f{t,xit),u{t))-x'{t))]dt. (45) 
Jo 

The optimal solutions for the optimal control problem can be found from the Euler-Lagrange 
equations for /. 

Theorem 5.1 (Euler-Lagrange Equation) If x is an optimal solution of J^F{t,x{t),u{t)) dt, then 
X must also be a solution of the problem (E-L) 

div{F^{t,x{t),Vx{t))) = F^^{t,x{t),Vx{t)) £ n, x = xq £ dn, (46) 

for the variables (x,'SJx) with either prescribed or transversality conditions applied on dO,. 
The Euler-Lagrange equations for the optimal control problem can be found first by defining 

G{t,u,p,x,u',p,x') = F{t,x,u) + p- {f{t,x,u) -x'), (47) 



where we note that the right-hand side of (47) has no explicit dependence on either u' or p'. 
Then the system can be written as 

d^dG__dG_ d^dG_dG d^dG _ dG 
dt dx' dx ' dt du' du ' dt dp' dp ' 

and, therefore, 



(48) 



, dF \ df . . /, X 

-P = -^(t, X, u) + P-r^{t, X, u), (49a) 

dF f 

= — (x, u, t) + p^{t, X, u), (49b) 

= x' - f{x,u,t). (49c) 

Defining the Hamiltonian of the system as H = F + pf, we can write 

, dH , , 

- = - & ■ f™") 

H{u) = mm H{v), (50b) 

X = /(t, X, u). (50c) 
where K is the set of all admissible controls. 



This provides a system of first-order differential equations, (50a) and (50c), for which we need 
boundary conditions, while (50b) is an algebraic constraint, which includes the end points in the 
case of a finite domain. We have at least one boundary condition for the equation involving x', while 
the boundary condition for is completed with the transversality condition or natural boundary 
condition. 

Theorem 5.2 (Transversality Condition) If at a given endpoint (initial or final) we have a con- 
dition on the state, we do not enforce the corresponding transversality condition, but if the state is 
free, then the transversality condition p = at the given endpoint must be taken into account J[^. 

This section provides a framework for an optimal control problem involving the diffusion of chemical 
species across an electrode domain. 
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5.3 Single Species Optimization 

In order to build an optimization routine without the use of built-in functions, a simpler model 
was introduced that contained only one differential equation that represents a single species dif- 
fusing across the domain. Instead of minimizing the current, the objective function consisted of 
maximizing the reaction rate balanced with a function representing the costs and durability of the 
electrode. This type of optimization would ensure longevity for the cathode material as well as 
maximizing the current produced. 

Fick's Laws were again used to produce the reaction-diffusion modeQ 

-D^ = h{c) = -ac, (51) 

where D is the control parameter exercised on the system. This parameter will vary across the 
domain but remains constant to first order and we can avoid taking the derivative with respect 
to the position in order to complete the optimization. The function h{c) represents the reaction 
rate and a contains other variables of the system. Since we are using only one equation, we can 
represent the reactivity, a, as a single constant. 



The function g{D), shown in Figure 15 for different values of a, represents the costs and 
durability of the cathode material 

g{D) = a{D-Dof. (52) 

The objective function will measure how good the control is based upon the criteria discussed 
above. It is given by 

I(c, D)= [ Xh{c) + (1 - \)g{D) dx. (53) 



Jo 

Here, A can be between and 1, but chosen to be 0.5. If A is 1, it turns out that the objective 
function cannot be minimized as the optimization fails. If A is 0, the objective function is no longer 
maximizing the reaction rate. 

The parameter a is increased to prevent the diffusion coefficient from becoming too large or 
too small in the optimization process. If we have a large diffusion coefficient, the porous medium 
no longer has an effect on the diffusion across the domain. It also means that the electrode is 
more porous and less stable mechanically. Note that D cannot exceed the diffusivity in bulk gas. 
A small diffusion coefficient would require smaller pores which are difficult to manufacture for the 
same volume fraction of species, resulting in high costs. Therefore, as the optimization routine 
processes, the diffusion coefficient remains in the vicinity of the minimum value, Dq. Here, we 
choose Do = 1- 

The boundary conditions are both considered at the origin, which represents the electrode/elec- 
trolyte interface where the gas cannot flow across 

dc 

c(0) = Co, —(0) = 0. (54) 



As the species diffuses from right to left in the domain (Figure 16), we expect the concentration 



to reach a specified value, Cq, while maintaining no flux at the left boundary into the electrolyte. 

^Strictly speaking, we need to solve — ^ (^(*)^) ~ ~ lE' 1^ ~ ~ ^^c. Hence, we shall see that for large a 

in (52 1, the term can be neglected. As a decreases, the term gains in importance. Future work will include 



the optimization of the full model — ^ (-^(^)5i) ~ possibly by utilizing the routine fmincon in MATLAB or 



other commercial software. 
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Figure 15: Function to be minimized representing the costs and durability of the cathode materials. If the 
diffusivity, D, is increased, the porous medium becomes less durable. If the diffusivity is decreased, the pore 
size becomes much smaller and is harder to manufacture for the same volume fraction of species. In order 
to stay near the optimal value (minimum D), the parameter a is increased. 



The channel is considered to be at the right boundary. The Dirichlet condition (c(0) = Cq), which 
is typically found at the channel (x = L) is moved to the electrolyte interface so as to have a 
well-defined optimization problem. 

Using the method described in Section 5.2 we start by formulating the second-order differential 
equation, Eq. (51 ), as two first-order differential equations with corresponding boundary conditions, 

a;i(0) = co, (55) 
X2{0) = 0, (56) 



4 



X2 



X2, 

D 



where the Hamiltonian is given by. 



H = Xh{xi) + (1 - X)g{D) + px2 



qhjxi) 

D ' 



(57) 



Here, p and q are co-states of x[ and Xg. 

Using the Hamiltonian, a system of equations is defined based upon the Euler-Lagrange system. 
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Figure 16: Domain schematic for the simpHfied model (51)-(54l 



Eqs. (50a)-(50c 



P 



dH 

dxi 
dH 

dx2 



Xh'{xi) 



qh'jxi) 
D 



-h'ixi) (A 



-P, 



H{D) = m.\nH{v) = min 



D\D - Do) 



Xh{xi) + (1 - X)g{v) +px2 
qh{xi) 



D 



qh{xi) 

V 



a 



(58a) 
(58b) 
(58c) 
(58d) 



where we used A = 1/2. 

The transversality conditions for this problem are given at the channel, 

p{L) = q{L) = 0. 



(59) 



The solution will generally require numerical tools. However, the problem will first be simplified 
in order to achieve an analytical solution of the system of four ODEs, including the algebraic 



equation (58d). 



In our simple model, lei a = cq = L = Dq = 1 and A 



Then Eqs. (55), (56), (58a) and 



(58b) become 



4 



P 



X2, 
Xl 

D' 
i_ 

2 D' 



Q = -P, 
with boundary conditions 

xi(0) = l, X2(0) = 0, 



p{l) 



0. 



(60a) 
(60b) 

(60c) 
(60d) 

(61) 



In Eqs. (60b) and (60c), D is determined by Eq. ( |58d| ), which couples all four equations (60a)- 
(60d). We see that as q ^ cxd, D — > 1 across the domain since D = is not permissible. Setting 



D = 1 decouples the four equations into 2 pairs of equations. This is acceptable as a first-order 
approximation to find an analytical solution since we will see that this is consistent for small a as 
well. Therefore, we first solve 



/ 



X2, 
Xl, 



XliO) = 1, 

X2(0) = 0, 



(62a) 
(62b) 
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Figure 17: The diffusion coefficient will always be a positive- valued number, which is determined by 
Eq. (58d). The smallest value of D will occur at the electrode/electrolyte boundary where the gas can 



no longer diffuse. The minimum value, D, is found where there exists a local minimum in Eq. (58dl as 
shown on the graph. 



which yields the general solutions 

xi{x) = Ci sinh(x) + C2 cosh(x), 
X2ix) = Ci cosh(x) + C2 sinh(x). 

Applying the boundary conditions gives the final solutions 

xi{x) = cosh(x), 
X2{x) = sinh(a;). 



Similarly, we solve the remaining two ODEs 



P 



q = -p, 

This gives the following solutions 

p{x) = Ci sinh(x) + C2 cosh(x), 

q{x) = - — (Ci cosh(x) + C2 sinh(x)). 



p(l) = 0, 
q{l)=0. 



(63a) 
(63b) 



(64a) 
(64b) 



(65a) 
(65b) 

(66a) 
(66b) 
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Figure 18: As the value of a decreases in g{D), the 
solution begins to diverge as the local minimum of 
D is passed. 



Figure 19: As the value of ac is approached, the 
CPU time begins to increase exponentially as it 
becomes much more difficult to reach a steady state 
solution numerically. 



Applying the boundary conditions leads to the two conditions 
p(l) = = Ci sinh(l) + C2 cosh(l), 
q{l) = = ^ - Ci cosh(l) - C2 sinh(l) 



(67a) 
(67b) 



with solutions 2Ci = cosh(l) and 2C2 
that 



-sinh(l). Applying these constants to p and q, one finds 



1 1 
p{x) = - cosh(l) sinh(x) — - sinh(l) cosh(x), 

q{x) ~ 2 ~ 2 cosh(x) + - sinh(l) sinh(x). 



(68a) 
(68b) 



The solutions to the system of ODEs in the limit a — > 00 are given in equations (64a)-(64b) and 



(68a)-(68b). 



In order to solve this problem numerically for general a, a shooting method is used. Based 
upon the analytical results found here for a — > 00, an initial guess can be obtained from p(0) and 



p{0) 



sinh(l) < 0, 



g(0) = -(1 -cosh(l)) < 0. 



(69a) 



Hence, we start in the 3rd quadrant of the (p, q) — plane and need to end up at the origin when 
x = l. 

A secant method is used to update the guess for the next iteration for a which needs two initial 
guesses. The second guess is chosen to be within 1% of the first initial guess so as not to stray away 
from the analytical results. By choosing the solution of the previous a as the next initial guess, a 
can be varied and decreased step-by-step, and the solutions can be found. Note that we are strictly 
shooting in four dimensions, equivalent to a four- dimensional surface embedded in five dimensions 



owing to the algebraic constraint ( 58d ) . The algebraic equation is solved using a built-in MATLAB 



root solver, fzero, after every iteration. 
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Figure 20: The concentration profile across the domain for decreasing values of a. As the value of is 
approached, the solution remains within the same order of magnitude. 



In order for the solution to converge, the right-hand side of Eq. (58d) cannot drop below the 



local minimum of the left-hand side or else D would be negative which is unphysical. The minimum, 
D, is found from. 



d 

dD 



D\D-Do) 



3D^ 



2DDq 
D 



0, 

0, 
2 



-Do, 



(70a) 
(70b) 
(70c) 



where Dq = 1 is chosen so as to simplify the system. The value, Dq, will occur at the right boundary 
since we have q{l) = in (58d). The diffusivity as a function of x is found using a root finder at 
each point, using (58d) and the solutions for q and xi. The minimum of il?o will occur close to 



some critical value of a, which is the lowest possible value where the numerical solution does not 
exist, and represents the smallest possible value of D and a, as seen in Figure 17 The minimum 



value of D decreases monotonically with a. The difference between the numerical minimum value, 
min{D}, and the analytical value, (|) (| — l) = — ^, is shown in Figure 
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Numerically to four decimal places, the critical value of a was computed as 4.1558. As the 



critical value is approached from above, the CPU time increases without bound as in Figure [19 
As a — > Oc, the solution approaches the graph shown in Figure 20 



Finally, we would like to gauge whether dropping the term 



The presence of this term would modify Eq. (51) to give 



in deriving (51) is justified. 



dD dc 

dx dx 



D 



d^ 
dx"^ 



-ac. 



(71) 
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Figure 21: Comparison of the two terms on the left-hand side of (71 1. For large a, the ratio is large and 
the model (51 1 seems justified. As a decreases, the extra term in the full model (71) becomes increasingly 
important. 



Dc" 
D'c' 



However, the 



For a = 1, the two terms on the left-hand side can be compared by their ratio 
presence of this additional term renders our optimization method non-applicable. 

Regardless of this fact, we are plotting in Figure [2T] this ratio for different values of a, approx- 
imated by our solution of the optimization problem xi{x), X2{x), D{x) and so 



(72) 



Dc" 




Dx'^ 


D'd 




D'X2 



We see that for large a, our original simplification in deriving (51) seems justified. However, as 
ttc we can observe that the term plays an increasingly important role and, in principle. 



a 



the full model needs to be optimized, i.e., Eq. (71) 



6 Conclusion and Future Work 

Three different models of diffusion across the cathode of a MCFC based on Fickian diffusion, 
convection-diffusion and multicomponent diffusion with convection were studied. It has been shown 
that the results can differ significantly, depending on system parameters. 

The convection-diffusion model shows that for standard values found from the literature, the 
convective flux dominates the gas species' total flux across the domain. This is a significant result 
since this phenomenon has not been taken into account by other researchers such as White et al. |2] 
and should be included in any model of mass transport for a MCFC. Since the convective flux is 
the dominant term, an approximation using only Fickian diffusion with an effective diffusivity is 
not possible as diffusion will only model a decrease in concentration towards the electrolyte. 
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The Maxwell-Stefan equations for multicomponent diffusion take into consideration the mo- 
mentum losses to due the interactions between particles of different species. Combining it with the 
convective flux gives a more detailed view of the mass transport inside the cathode. For the stan- 
dard values from the literature, the results of the multicomponent model follow those of the simpler 
convection-diffusion model with very high accuracy. In other words, the convection-diffusion model 
is a good approximation of the more complex binary diffusion model. At low values of the per- 
meability, the convective flux approaches zero, which allows the differences between the models 
to become apparent since the diffusive fluxes do (simple diffusion) or do not (multi-component 
diffusion) add up to zero. 

While performing the numerical simulations for values of the liquid conductivity obtained from 
the literature, a steady-state solution could not be achieved for small values. By using a simpler 
problem, an analytical solution was derived showing that there exists a critical value for the liquid 
conductivity below which steady-state solutions do not exist. This is another significant result as it 
shows that not only convection needs to be included in the White et al. [2] model but the parameter 
values are inconsistent. 

The optimization of the electrode for a fuel cell is of great importance as it will aid in the 
manufacturing of cell components to improve cell performance as well as increase life time. An 
analytical solution has been derived for a simplified system that involves the optimal profile of the 
diffusivity across the domain as a control on the cathode. The diffusivity can be varied in the more 
comprehensive models using the porosity profile across the domain. Using the data obtained from 
the optimization routine, both analytical and numerical, the porosity may be manufactured so as 
to improve the performance and life time of the cell. 

Future work depends upon the availability of parameters that will help to better understand the 
convection-diffusion results. As the fuel cell begins to reach market acceptance, the optimization 
of the fuel cell will grow in importance and will aid in the manufacturing of cell components to 
further increase the viability of the MCFC. A more realistic optimization of the cathode transport 
processes needs to be investigated. 
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